####### INTER-CASTE INTERACTIONS #############
options(scipen=999)
reds = reds %>%
  group_by(stateid, districtid, tehsilid, blockid, villageid) %>%
  mutate(vwreserv_gp1 = ifelse(vwreserv==1 & period_new==1, 1, 0) %>% max(na.rm = T)) %>%
  mutate(vwreserv_gp2 = ifelse(vwreserv==1 & period_new==2, 1, 0) %>% max(na.rm = T)) %>%
  mutate(vwreserv_gp3 = ifelse(vwreserv==1 & period_new==3, 1, 0) %>% max(na.rm = T)) %>%
  mutate(vstreserv_gp1 = ifelse(vstreserv==1 & period_new==1, 1, 0) %>% max(na.rm = T)) %>%
  mutate(vstreserv_gp2 = ifelse(vstreserv==1 & period_new==2, 1, 0) %>% max(na.rm = T)) %>%
  mutate(vstreserv_gp3 = ifelse(vstreserv==1 & period_new==3, 1, 0) %>% max(na.rm = T)) %>%
  mutate(vscreserv_gp1 = ifelse(vscreserv==1 & period_new==1, 1, 0) %>% max(na.rm = T)) %>%
  mutate(vscreserv_gp2 = ifelse(vscreserv==1 & period_new==2, 1, 0) %>% max(na.rm = T)) %>%
  mutate(vscreserv_gp3 = ifelse(vscreserv==1 & period_new==3, 1, 0) %>% max(na.rm = T)) %>%
  mutate(vwr_vstreserv_gp1 = ifelse(vwreserv==1 & vstreserv==1 & period_new==1, 1, 0) %>% max(na.rm = T)) %>%
  mutate(vwr_vstreserv_gp2 = ifelse(vwreserv==1 & vstreserv==1 & period_new==2, 1, 0) %>% max(na.rm = T)) %>%
  mutate(vwr_vstreserv_gp3 = ifelse(vwreserv==1 & vstreserv==1 & period_new==3, 1, 0) %>% max(na.rm = T)) %>%
  mutate(vwreserv_gp12 = ifelse(vwreserv_gp1==1 | vwreserv_gp2==1, 1, 0)) %>%
  mutate(vstreserv_gp12 = ifelse(vstreserv_gp1==1 | vstreserv_gp2==1, 1, 0)) %>%
  mutate(vscreserv_gp12 = ifelse(vscreserv_gp1==1 | vscreserv_gp2==1, 1, 0)) %>%
  mutate(vscstreserv_gp1 = ifelse(vscstreserv==1 & period_new==1, 1, 0) %>% max(na.rm = T)) %>%
  mutate(vscstreserv_gp2 = ifelse(vscstreserv==1 & period_new==2, 1, 0) %>% max(na.rm = T)) %>%
  mutate(vscstreserv_gp3 = ifelse(vscstreserv==1 & period_new==3, 1, 0) %>% max(na.rm = T)) %>%
  mutate(vscstreserv_gp12 = ifelse(vscreserv_gp12==1 | vstreserv_gp12==1, 1, 0)) %>%
  mutate(vwrXvst_gp12 = ifelse(vwrXvst_gp1==1 | vwrXvst_gp2==1, 1, 0)) %>%
  
  mutate(vwreserv_prior = case_when(
    (vwreserv_gp1==1|vwreserv_gp2==1) & period_new==3 ~ 1,
    (vwreserv_gp1==0&vwreserv_gp2==0) & period_new==3 ~ 0,
    (vwreserv_gp1==1) & period_new==2 ~ 1,
    (vwreserv_gp1==0) & period_new==2 ~ 0,
    period_new==1 ~ 0)) %>%
  mutate(vstreserv_prior = case_when(
    (vstreserv_gp1==1|vstreserv_gp2==1) & period_new==3 ~ 1,
    (vstreserv_gp1==0&vstreserv_gp2==0) & period_new==3 ~ 0,
    (vstreserv_gp1==1) & period_new==2 ~ 1,
    (vstreserv_gp1==0) & period_new==2 ~ 0,
    period_new==1 ~ 0)) %>%
  mutate(vscreserv_prior = case_when(
    (vscreserv_gp1==1|vscreserv_gp2==1) & period_new==3 ~ 1,
    (vscreserv_gp1==0&vscreserv_gp2==0) & period_new==3 ~ 0,
    (vscreserv_gp1==1) & period_new==2 ~ 1,
    (vscreserv_gp1==0) & period_new==2 ~ 0,
    period_new==1 ~ 0)) %>% 
  mutate(vscstreserv_prior = case_when(
    (vscstreserv_gp1==1|vscstreserv_gp2==1) & period_new==3 ~ 1,
    (vscstreserv_gp1==0&vscstreserv_gp2==0) & period_new==3 ~ 0,
    (vscstreserv_gp1==1) & period_new==2 ~ 1,
    (vscstreserv_gp1==0) & period_new==2 ~ 0,
    period_new==1 ~ 0)) 

# Matched datasets
reds.st.matched = reds.st.matched %>%
  group_by(stateid, districtid, tehsilid, blockid, villageid) %>%
  mutate(vwreserv_gp1 = ifelse(vwreserv==1 & period_new==1, 1, 0) %>% max(na.rm = T)) %>%
  mutate(vwreserv_gp2 = ifelse(vwreserv==1 & period_new==2, 1, 0) %>% max(na.rm = T)) %>%
  mutate(vwreserv_gp3 = ifelse(vwreserv==1 & period_new==3, 1, 0) %>% max(na.rm = T)) %>%
  mutate(vstreserv_gp1 = ifelse(vstreserv==1 & period_new==1, 1, 0) %>% max(na.rm = T)) %>%
  mutate(vstreserv_gp2 = ifelse(vstreserv==1 & period_new==2, 1, 0) %>% max(na.rm = T)) %>%
  mutate(vstreserv_gp3 = ifelse(vstreserv==1 & period_new==3, 1, 0) %>% max(na.rm = T)) %>%
  mutate(vwr_vstreserv_gp1 = ifelse(vwreserv==1 & vstreserv==1 & period_new==1, 1, 0) %>% max(na.rm = T)) %>%
  mutate(vwr_vstreserv_gp2 = ifelse(vwreserv==1 & vstreserv==1 & period_new==2, 1, 0) %>% max(na.rm = T)) %>%
  mutate(vwr_vstreserv_gp3 = ifelse(vwreserv==1 & vstreserv==1 & period_new==3, 1, 0) %>% max(na.rm = T)) %>%
  mutate(vwreserv_gp12 = ifelse(vwreserv_gp1==1 | vwreserv_gp2==1, 1, 0)) %>%
  mutate(vstreserv_gp12 = ifelse(vstreserv_gp1==1 | vstreserv_gp2==1, 1, 0)) %>%
  mutate(vwreserv_prior = case_when(
    (vwreserv_gp1==1|vwreserv_gp2==1) & period_new==3 ~ 1,
    (vwreserv_gp1==0&vwreserv_gp2==0) & period_new==3 ~ 0,
    (vwreserv_gp1==1) & period_new==2 ~ 1,
    (vwreserv_gp1==0) & period_new==2 ~ 0,
    period_new==1 ~ 0)) %>%
  mutate(vstreserv_prior = case_when(
    (vstreserv_gp1==1|vstreserv_gp2==1) & period_new==3 ~ 1,
    (vstreserv_gp1==0&vstreserv_gp2==0) & period_new==3 ~ 0,
    (vstreserv_gp1==1) & period_new==2 ~ 1,
    (vstreserv_gp1==0) & period_new==2 ~ 0,
    period_new==1 ~ 0)) 

reds.scst.matched = reds.scst.matched %>%
  group_by(stateid, districtid, tehsilid, blockid, villageid) %>%
  mutate(vwreserv_gp1 = ifelse(vwreserv==1 & period_new==1, 1, 0) %>% max(na.rm = T)) %>%
  mutate(vwreserv_gp2 = ifelse(vwreserv==1 & period_new==2, 1, 0) %>% max(na.rm = T)) %>%
  mutate(vwreserv_gp3 = ifelse(vwreserv==1 & period_new==3, 1, 0) %>% max(na.rm = T)) %>%
  mutate(vscstreserv_gp1 = ifelse(vscstreserv==1 & period_new==1, 1, 0) %>% max(na.rm = T)) %>%
  mutate(vscstreserv_gp2 = ifelse(vscstreserv==1 & period_new==2, 1, 0) %>% max(na.rm = T)) %>%
  mutate(vscstreserv_gp3 = ifelse(vscstreserv==1 & period_new==3, 1, 0) %>% max(na.rm = T)) %>%
  mutate(vwreserv_none = ifelse(vwreserv_gp1==0 | vwreserv_gp2==0, 1, 0)) %>%
  mutate(vwreserv_gp12 = ifelse(vwreserv_gp1==1 | vwreserv_gp2==1, 1, 0)) %>%
  mutate(vscstreserv_gp12 = ifelse(vscstreserv_gp1==1 | vscstreserv_gp2==1, 1, 0)) %>%
  mutate(vwreserv_prior = case_when(
    (vwreserv_gp1==1|vwreserv_gp2==1) & period_new==3 ~ 1,
    (vwreserv_gp1==0&vwreserv_gp2==0) & period_new==3 ~ 0,
    (vwreserv_gp1==1) & period_new==2 ~ 1,
    (vwreserv_gp1==0) & period_new==2 ~ 0,
    period_new==1 ~ 0)) %>%
  mutate(vscstreserv_prior = case_when(
    (vscstreserv_gp1==1|vscstreserv_gp2==1) & period_new==3 ~ 1,
    (vscstreserv_gp1==0&vscstreserv_gp2==0) & period_new==3 ~ 0,
    (vscstreserv_gp1==1) & period_new==2 ~ 1,
    (vscstreserv_gp1==0) & period_new==2 ~ 0,
    period_new==1 ~ 0)) 


# Set up samples
reds_p3 = reds %>% filter(period_new==3) # Data only from period 3
reds_p3_rand = reds_p3 %>% filter(nr_short!=1) # Data only from states with random assignment, period 3



#################### TABLE D.5 #############
reds_p3 = reds_p3 %>%
  mutate(woman_pradhanXST_prad = case_when(
    woman_pradhan==1 & ST_pradhan==1 ~ 1, 
    woman_pradhan==0 | ST_pradhan==0 ~ 0, 
    TRUE ~ NA_real_
  )) %>%
  mutate(vwreservXvstreserv = case_when(
    vwreserv==1 & vstreserv==1 ~ 1, 
    vwreserv==0 | vstreserv==0 ~ 0, 
    TRUE ~ NA_real_
  ))

controls = "top20land + birthyr + fem + Caste_ST + Caste_SC + vwreserv_gp12 + vstreserv_gp12"

fs_st = felm(paste("caste_interact_bi_comp ~ prop_STnow + ", controls, "| tehsilid |  (woman_pradhan|ST_pradhan|I(woman_pradhan*ST_pradhan) ~ vwreserv+vstreserv+I(vwreserv*vstreserv)) | villageid") %>% as.formula,
             data = reds_p3_rand)
fs_scst = felm(paste("caste_interact_bi_comp ~ prop_STnow + ", controls, "| tehsilid |  (woman_pradhan|ST_SC_pradhan|I(woman_pradhan*ST_SC_pradhan) ~ vwreserv+vscstreserv+I(vwreserv*vscstreserv)) | villageid") %>% as.formula,
               data = reds_p3_rand)

summary(fs_st$stage1, lhs = "woman_pradhan") 
summary(fs_scst$stage1, lhs = "woman_pradhan")
summary(fs_st$stage1, lhs = "ST_pradhan")
summary(fs_scst$stage1, lhs = "ST_SC_pradhan")
summary(fs_st$stage1, lhs = "I(woman_pradhan * ST_pradhan)")
summary(fs_scst$stage1, lhs = "I(woman_pradhan * ST_SC_pradhan)")

row1 <- c(tidy(fs_st$stage1)[9,"estimate"], tidy(fs_scst$stage1)[9,"estimate"], 
          tidy(fs_st$stage1)[20,"estimate"],tidy(fs_scst$stage1)[20,"estimate"], 
          tidy(fs_st$stage1)[31,"estimate"],tidy(fs_scst$stage1)[31,"estimate"])
row2 <- c(tidy(fs_st$stage1)[9,"std.error"], tidy(fs_scst$stage1)[9,"std.error"], 
          tidy(fs_st$stage1)[20,"std.error"],tidy(fs_scst$stage1)[20,"std.error"], 
          tidy(fs_st$stage1)[31,"std.error"],tidy(fs_scst$stage1)[31,"std.error"])
row3 <- c(tidy(fs_st$stage1)[10,"estimate"], tidy(fs_scst$stage1)[10,"estimate"], 
          tidy(fs_st$stage1)[21,"estimate"],tidy(fs_scst$stage1)[21,"estimate"], 
          tidy(fs_st$stage1)[32,"estimate"],tidy(fs_scst$stage1)[32,"estimate"])
row4 <- c(tidy(fs_st$stage1)[10,"std.error"], tidy(fs_scst$stage1)[10,"std.error"], 
          tidy(fs_st$stage1)[21,"std.error"],tidy(fs_scst$stage1)[21,"std.error"], 
          tidy(fs_st$stage1)[32,"std.error"],tidy(fs_scst$stage1)[32,"std.error"])
row5 <- c(tidy(fs_st$stage1)[11,"estimate"], tidy(fs_scst$stage1)[11,"estimate"], 
          tidy(fs_st$stage1)[22,"estimate"],tidy(fs_scst$stage1)[22,"estimate"], 
          tidy(fs_st$stage1)[33,"estimate"],tidy(fs_scst$stage1)[33,"estimate"])
row6 <- c(tidy(fs_st$stage1)[11,"std.error"], tidy(fs_scst$stage1)[11,"std.error"], 
          tidy(fs_st$stage1)[22,"std.error"],tidy(fs_scst$stage1)[22,"std.error"], 
          tidy(fs_st$stage1)[33,"std.error"],tidy(fs_scst$stage1)[33,"std.error"])

header <- " & Woman Pradhan & Woman Pradhan & ST Pradhan & ST/SC Pradhan & Woman ST Pradhan & Woman ST/SC Pradhan \\\\"
row1_text <- sprintf("W Quota & %.3f$^{***}$ & %.3f$^{***}$ & %.3f$^{**}$ & %.3f$^{***}$ & %.3f$^{***}$ & %.3f$^{***}$ \\\\",
                     row1[1], row1[2], row1[3], row1[4], row1[5], row1[6])
row2_text <- sprintf(" & (%.3f) & (%.3f) & (%.3f) & (%.3f) & (%.3f) & (%.3f) \\\\",
                     row2[1], row2[2], row2[3], row2[4], row2[5], row2[6])
row3_text <- sprintf("Caste Quota & %.3f$^{***}$ & %.3f & %.3f$^{***}$ & %.3f$^{***}$ & %.3f$^{**}$ & %.3f \\\\",
                     row3[1], row3[2], row3[3], row3[4], row3[5], row3[6])
row4_text <- sprintf(" & (%.3f) & (%.3f) & (%.3f) & (%.3f) & (%.3f) & (%.3f) \\\\",
                     row4[1], row4[2], row4[3], row4[4], row4[5], row4[6])
row5_text <- sprintf("W X C Quota & %.3f & %.3f & %.3f & %.3f$^{*}$ & %.3f$^{***}$ & %.3f$^{**}$ \\\\",
                     row5[1], row5[2], row5[3], row5[4], row5[5], row5[6])
row6_text <- sprintf(" & (%.3f) & (%.3f) & (%.3f) & (%.3f) & (%.3f) & (%.3f) \\\\",
                     row6[1], row6[2], row6[3], row6[4], row6[5], row6[6])

# Combine all rows into a full LaTeX table
# Open a sink connection to the file
sink(paste0(tab.out, "tableD5.tex"))

# Print the LaTeX table using cat()
cat("\\begin{table}[ht]\n\\centering\n")
cat("\\begin{tabular}{lcccccc}\n")
cat(header, "\n")
cat("\\hline\n")
cat(row1_text, "\n")
cat(row2_text, "\n")
cat(row3_text, "\n")
cat(row4_text, "\n")
cat(row5_text, "\n")
cat(row6_text, "\n")
cat("\\hline\n")
cat("Caste Quota Type & ST & ST/SC & ST & ST/SC & ST & ST/SC\n")
cat("Observations & 83329 & 83329 & 83329 & 83329 & 83329 & 83329\n")
cat("Sub-district FE & Yes & Yes & Yes & Yes & Yes & Yes\n")
cat("Controls & Yes & Yes & Yes & Yes & Yes & Yes\n")
cat("\\hline\n")
cat("\\end{tabular}\n")
cat("\\caption{Your Caption Here}\n")
cat("\\label{tab:yourlabel}\n")
cat("\\end{table}\n")

# Close the sink connection so that output goes back to the console
sink()



###### TABLE 2 ######

base = "caste_interact_bi_comp ~ vwreserv*vstreserv"
controls = "top20land + birthyr + fem + Caste_ST + Caste_SC + vwreserv_gp12 + vstreserv_gp12"

reds.lm3 = paste(base, "+", controls, "+", "prop_STnow") %>% as.formula() %>%
  lm_robust(data = reds_p3_rand,
            fixed_effects = ~districtid, clusters = villageid, se_type = "stata") 

reds.lm4 = paste(base, "+", controls, "+", "prop_STnow") %>% as.formula() %>%
  lm_robust(data = reds_p3_rand,
            fixed_effects = ~tehsilid, clusters = villageid, se_type = "stata") 

reds.lm5 = paste(base, "+", controls, "+ prop_STnow") %>% as.formula() %>%
  lm_robust(data = reds.st.matched %>% filter(period_new==3 & nr_short!=1),
            fixed_effects = ~districtid,
            clusters = villageid, se_type = "stata") 

base = "caste_interact_bi_comp ~ vwreserv*vscstreserv"
controls = "top20land + birthyr + fem + Caste_ST + Caste_SC + vwreserv_gp12 + vscstreserv_gp12"

reds.lm6 = paste(base, "+", controls, "+", "prop_SCnow + prop_STnow") %>% as.formula() %>%
  lm_robust(data = reds_p3_rand,
            fixed_effects = ~tehsilid, clusters = villageid, se_type = "stata") 

controls = "top20land + birthyr + fem + Caste_ST + Caste_SC + vwreserv_gp12 + vscstreserv_gp12"
reds.lm7 = paste(base, "+", controls, "+", "prop_SCnow + prop_STnow") %>% as.formula() %>%
  lm_robust(data = reds.scst.matched %>% filter(period_new==3 & nr_short!=1),
            fixed_effects = ~districtid, clusters = villageid, se_type = "stata")


### Analyzing ST and ST/SC women
base = "caste_interact_bi_comp ~ vwreserv*vstreserv"
controls = "top20land + birthyr + vwreserv_gp12 + vstreserv_gp12"
reds.lm8 =  paste(base, "+", controls, "+ prop_STnow") %>% as.formula() %>%
  lm_robust(data = reds.st.matched %>% filter(period_new==3 & nr_short!=1 & fem==1 & Caste_ST==1),
            fixed_effects = ~districtid,
            clusters = villageid, se_type = "stata") 

base = "caste_interact_bi_comp ~ vwreserv*vscstreserv"
controls = "top20land + birthyr + vwreserv_gp12 + vscstreserv_gp12"
reds.lm9 = paste(base, "+", controls, "+", "prop_SCnow + prop_STnow") %>% as.formula() %>%
  lm_robust(data = reds.scst.matched %>% filter(period_new==3 & nr_short!=1 & fem==1 & (Caste_ST==1|Caste_SC==1)),
            fixed_effects = ~districtid, clusters = villageid, se_type = "stata")


means = c(
  mean(reds_p3_rand$caste_interact_bi_comp, na.rm = T),
  mean(reds_p3_rand$caste_interact_bi_comp, na.rm = T),
  mean(reds.st.matched$caste_interact_bi_comp[reds.st.matched$period_new==3 & reds.st.matched$nr_short!=1], na.rm = T),
  mean(reds_p3_rand$caste_interact_bi_comp, na.rm = T),
  mean(reds.scst.matched$caste_interact_bi_comp[reds.scst.matched$period_new==3 & reds.scst.matched$nr_short!=1], na.rm = T),
  mean(reds.st.matched$caste_interact_bi_comp[reds.st.matched$period_new==3 & reds.st.matched$nr_short!=1 & reds.st.matched$fem==1 & reds.st.matched$Caste_ST==1], na.rm = T),
  mean(reds.scst.matched$caste_interact_bi_comp[reds.scst.matched$period_new==3 & reds.scst.matched$nr_short!=1 & reds.scst.matched$fem==1 & (reds.scst.matched$Caste_ST==1|reds.scst.matched$Caste_SC==1)], na.rm = T)
)


screenreg(l = list(reds.lm3, reds.lm4, reds.lm5, 
                reds.lm6, reds.lm7, reds.lm8, 
                reds.lm9), include.ci = F, include.rmse = F,
       include.nobs = T,include.adjrs = FALSE, include.rsquared = F,
       digits = 3, include.nclust = F,
       omit.coef = c("top20land|birthyr|fem|Caste_ST|Caste_SC|vwreserv_gp12|vstreserv_gp12|prop_STnow|vscstreserv_gp12|prop_SCnow"),
       custom.coef.names = c("Women's Quota", "Caste Quota (ST)", "W X C (ST) Quota", 
                             "Caste Quota (ST/SC)", "W X Caste (ST/SC) Quota"),
       stars = c(0.01, 0.05, 0.1),
       custom.model.names = c("All", "All", "All", "All", "All", "ST Women", "ST/SC Women"),
       custom.gof.rows = list(
         "DV Mean" = means,
         "Caste Quota Type" = c("ST", "ST", "ST", "ST/SC", "ST/SC", "ST", "ST/SC"),
         "District FE" = c("Yes", "No", "Yes", "No", "Yes", "Yes", "Yes"),
         "Subdistrict FE" = c("No", "Yes", "No", "Yes", "No", "No", "No"),
         "Matched Data" = c("No", "No", "Yes", "No", "Yes", "Yes", "Yes"),
         "Controls" = c("Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes")), 
       file = paste0(tab.out, "table2.tex"))






###### TABLE E6, TABLE E7 ######

base = "caste_interact_bi_comp ~ vwreserv*vstreserv"
controls = "top20land + birthyr + fem + Caste_ST + Caste_SC + vwreserv_gp12 + vstreserv_gp12"

reds.lm1 = paste(base) %>% as.formula() %>%
  lm_robust(data = reds_p3, clusters = villageid, se_type = "stata") 

reds.lm2 = paste(base) %>% as.formula() %>%
  lm_robust(data = reds_p3_rand,
            fixed_effects = ~districtid, clusters = villageid, se_type = "stata") 

reds.lm3 = paste(base, "+", controls, "+", "prop_STnow") %>% as.formula() %>%
  lm_robust(data = reds_p3_rand,
            fixed_effects = ~districtid, clusters = villageid, se_type = "stata") 

reds.lm4 = paste(base, "+", controls, "+", "prop_STnow") %>% as.formula() %>%
  lm_robust(data = reds_p3_rand,
            fixed_effects = ~tehsilid, clusters = villageid, se_type = "stata") 



reds.lm5 = paste(base, "+", controls, "+ prop_STnow") %>% as.formula() %>%
  lm_robust(data = reds.st.matched %>% filter(period_new==3 & nr_short!=1),
            fixed_effects = ~districtid,
            clusters = villageid, se_type = "stata") 

base = "caste_interact_bi_comp ~ vwreserv*vscstreserv"
controls = "top20land + birthyr + fem + Caste_ST + Caste_SC + vwreserv_gp12 + vscstreserv_gp12"

reds.lm6 = paste(base, "+", controls, "+", "prop_SCnow + prop_STnow") %>% as.formula() %>%
  lm_robust(data = reds_p3_rand,
            fixed_effects = ~tehsilid, clusters = villageid, se_type = "stata") 

controls = "top20land + birthyr + fem + Caste_ST + Caste_SC + vwreserv_gp12 + vscstreserv_gp12"
reds.lm7 = paste(base, "+", controls, "+", "prop_SCnow + prop_STnow") %>% as.formula() %>%
  lm_robust(data = reds.scst.matched %>% filter(period_new==3 & nr_short!=1),
            fixed_effects = ~districtid, clusters = villageid, se_type = "stata")


means = c(
  mean(reds_p3$caste_interact_bi_comp, na.rm = T),
  mean(reds_p3_rand$caste_interact_bi_comp, na.rm = T),
  mean(reds_p3_rand$caste_interact_bi_comp, na.rm = T),
  mean(reds_p3_rand$caste_interact_bi_comp, na.rm = T),
  mean(reds.st.matched$caste_interact_bi_comp[reds.st.matched$period_new==3 & reds.st.matched$nr_short!=1], na.rm = T),
  mean(reds_p3_rand$caste_interact_bi_comp, na.rm = T),
  mean(reds.scst.matched$caste_interact_bi_comp[reds.scst.matched$period_new==3 & reds.scst.matched$nr_short!=1], na.rm = T)
)




screenreg(l = list(reds.lm1, reds.lm2, reds.lm3, 
                reds.lm4, reds.lm5, reds.lm6, 
                reds.lm7), include.ci = F, include.rmse = F,
       include.nobs = T,include.adjrs = FALSE, include.rsquared = F,
       digits = 3, include.nclust = F,
       omit.coef = c("Intercept|top20land|birthyr|fem|Caste_ST|Caste_SC|vwreserv_gp12|vstreserv_gp12|prop_STnow|vscstreserv_gp12|prop_SCnow"),
       custom.coef.names = c("Women's Quota", "Caste Quota (ST)", "W X C (ST) Quota", 
                             "Caste Quota (ST/SC)", "W X Caste (ST/SC) Quota"),
       stars = c(0.01, 0.05, 0.1),
       custom.model.names = c("All", "All", "All", "All", "All", "All", "All"),
       custom.gof.rows = list(
         "DV Mean" = means,
         "Caste Quota Type" = c("ST", "ST", "ST", "ST", "ST", "ST/SC", "ST/SC"),
         "District FE" = c("No", "Yes", "Yes", "No", "Yes", "No", "Yes"),
         "Subdistrict FE" = c("No", "No", "No", "Yes", "No", "Yes", "No"),
         "Matched Data" = c("No", "No", "No", "No", "Yes", "No", "Yes"),
         "Controls" = c("No", "No", "Yes", "Yes", "Yes", "Yes", "Yes")), 
       file = paste0(tab.out, "tableE6.tex"))

means_short = c(
  mean(reds_p3_rand$caste_interact_bi_comp, na.rm = T),
  mean(reds.st.matched$caste_interact_bi_comp[reds.st.matched$period_new==3 & reds.st.matched$nr_short!=1], na.rm = T),
  mean(reds_p3_rand$caste_interact_bi_comp, na.rm = T),
  mean(reds.scst.matched$caste_interact_bi_comp[reds.scst.matched$period_new==3 & reds.scst.matched$nr_short!=1], na.rm = T)
)




screenreg(l = list(reds.lm4, reds.lm5, reds.lm6, 
                   reds.lm7), include.ci = F, include.rmse = F,
          include.nobs = T,include.adjrs = FALSE, include.rsquared = F,
          digits = 3, include.nclust = F,
          omit.coef = c("Intercept"),
          custom.coef.names = c("Women's Quota", "Caste Quota (ST)",
                                "Top 20\\% Landowner", "Birth Year", 
                                "Female", "Caste: ST", "Caste: SC", "Prior W Quota", 
                                "Prior Caste (ST) Quota", "Prop. ST", 
                                "W X C (ST) Quota", 
                                "Caste Quota (ST/SC)", "Prior Caste Quota (ST/SC)",
                                "Prop. SC",
                                "W X Caste (ST/SC) Quota"),
          stars = c(0.01, 0.05, 0.1),
          custom.model.names = c("All", "All", "All", "All"),
          custom.gof.rows = list(
            "DV Means" = means_short,
            "Sample" = c("Full", "Full", "Full", "Full"),
            "Caste Quota Type" = c("ST", "ST", "ST/SC", "ST/SC"),
            "District FE" = c("No", "Yes", "No", "Yes"),
            "Subdistrict FE" = c("Yes", "No", "Yes", "No"),
            "Matched Data" = c("No", "Yes", "No", "Yes"),
            "Controls" = c("Yes", "Yes", "Yes", "Yes")), 
          file = paste0(tab.out, "tableE7.tex"))


####### TABLE E.8 ##########
base = "caste_interact_bi_comp ~ vwreserv*vstreserv"
# Drop control for gender because data is restricted to women
controls = "top20land + birthyr + Caste_ST + Caste_SC + vwreserv_gp12 + vstreserv_gp12"

reds.lm1 = paste(base) %>% as.formula() %>%
  lm_robust(data = reds_p3 %>% filter(fem==1), clusters = villageid, se_type = "stata") 

reds.lm2 = paste(base) %>% as.formula() %>%
  lm_robust(data = reds_p3_rand %>% filter(fem==1),
            fixed_effects = ~districtid, clusters = villageid, se_type = "stata") 

reds.lm3 = paste(base, "+", controls, "+", "prop_STnow") %>% as.formula() %>%
  lm_robust(data = reds_p3_rand %>% filter(fem==1),
            fixed_effects = ~districtid, clusters = villageid, se_type = "stata") 

reds.lm4 = paste(base, "+", controls, "+", "prop_STnow") %>% as.formula() %>%
  lm_robust(data = reds_p3_rand %>% filter(fem==1),
            fixed_effects = ~tehsilid, clusters = villageid, se_type = "stata") 

reds.lm5 = paste(base, "+", controls, "+ prop_STnow") %>% as.formula() %>%
  lm_robust(data = reds.st.matched %>% filter(period_new==3 & nr_short!=1 & fem==1),
            fixed_effects = ~districtid,
            clusters = villageid, se_type = "stata") 

base = "caste_interact_bi_comp ~ vwreserv*vscstreserv"
reds.lm6 = paste(base, "+", controls, "+", "prop_SCnow + prop_STnow") %>% as.formula() %>%
  lm_robust(data = reds_p3_rand %>% filter(fem==1),
            fixed_effects = ~tehsilid, clusters = villageid, se_type = "stata") 



controls = "top20land + birthyr + Caste_ST + Caste_SC + vwreserv_gp12 + vscstreserv_gp12"
reds.lm7 = paste(base, "+", controls, "+", "prop_SCnow + prop_STnow") %>% as.formula() %>%
  lm_robust(data = reds.scst.matched %>% filter(period_new==3 & nr_short!=1 & fem==1),
            fixed_effects = ~districtid, clusters = villageid, se_type = "stata")

means = c(
  mean(reds_p3$caste_interact_bi_comp[reds_p3$fem==1], na.rm = T),
  mean(reds_p3_rand$caste_interact_bi_comp[reds_p3_rand$fem==1], na.rm = T),
  mean(reds_p3_rand$caste_interact_bi_comp[reds_p3_rand$fem==1], na.rm = T),
  mean(reds_p3_rand$caste_interact_bi_comp[reds_p3_rand$fem==1], na.rm = T),
  mean(reds.st.matched$caste_interact_bi_comp[reds.st.matched$period_new==3 & reds.st.matched$nr_short!=1 & reds.st.matched$fem==1], na.rm = T),
  mean(reds_p3_rand$caste_interact_bi_comp[reds_p3_rand$fem==1], na.rm = T),
  mean(reds.scst.matched$caste_interact_bi_comp[reds.scst.matched$period_new==3 & reds.scst.matched$nr_short!=1 & reds.scst.matched$fem==1], na.rm = T)
)


screenreg(l = list(reds.lm1, reds.lm2, reds.lm3, 
                   reds.lm4, reds.lm5, reds.lm6, 
                   reds.lm7), include.ci = F, include.rmse = F,
          include.nobs = T,include.adjrs = FALSE, include.rsquared = F,
          digits = 3, include.nclust = F,
          omit.coef = c("Intercept|top20land|birthyr|fem|Caste_ST|Caste_SC|vwreserv_gp12|vstreserv_gp12|prop_STnow|vscstreserv_gp12|prop_SCnow"),
          custom.coef.names = c("Women's Quota", "Caste Quota (ST)", "W X C (ST) Quota", 
                                "Caste Quota (ST/SC)", "W X Caste (ST/SC) Quota"),
          stars = c(0.01, 0.05, 0.1),
          custom.model.names = c("Women", "Women", "Women", "Women", "Women", "Women", "Women"),
          custom.gof.rows = list(
            "DV Means" = means,
            "Caste Quota Type" = c("ST", "ST", "ST", "ST", "ST", "ST/SC", "ST/SC"),
            "District FE" = c("No", "Yes", "Yes", "No", "Yes", "No", "Yes"),
            "Subdistrict FE" = c("No", "No", "No", "Yes", "No", "Yes", "No"),
            "Matched Data" = c("No", "No", "No", "No", "Yes", "No", "Yes"),
            "Controls" = c("No", "No", "Yes", "Yes", "Yes", "Yes", "Yes")), 
          file = paste0(tab.out, "tableE8.tex"))


### TABLE E.9 ######
base = "caste_interact_bi_comp ~ vwreserv*vstreserv"
controls = "top20land + birthyr + vwreserv_gp12 + vstreserv_gp12"

reds.lm1 = paste(base) %>% as.formula() %>%
  lm_robust(data = reds_p3 %>% filter(fem==1 & Caste_ST==1), clusters = villageid, se_type = "stata") 

reds.lm2 = paste(base) %>% as.formula() %>%
  lm_robust(data = reds_p3_rand %>% filter(fem==1 & Caste_ST==1),
            fixed_effects = ~districtid, clusters = villageid, se_type = "stata") 

reds.lm3 = paste(base, "+", controls, "+", "prop_STnow") %>% as.formula() %>%
  lm_robust(data = reds_p3_rand %>% filter(fem==1 & Caste_ST==1),
            fixed_effects = ~districtid, clusters = villageid, se_type = "stata") 

reds.lm4 = paste(base, "+", controls, "+", "prop_STnow") %>% as.formula() %>%
  lm_robust(data = reds_p3_rand %>% filter(fem==1 & Caste_ST==1),
            fixed_effects = ~tehsilid, clusters = villageid, se_type = "stata") 

reds.lm5 = paste(base, "+", controls, "+ prop_STnow") %>% as.formula() %>%
  lm_robust(data = reds.st.matched %>% filter(period_new==3 & nr_short!=1 & fem==1 & Caste_ST==1),
            fixed_effects = ~districtid,
            clusters = villageid, se_type = "stata") 

base = "caste_interact_bi_comp ~ vwreserv*vscstreserv"
reds.lm6 = paste(base, "+", controls, "+", "prop_SCnow + prop_STnow") %>% as.formula() %>%
  lm_robust(data = reds_p3_rand %>% filter(fem==1 & (Caste_ST==1|Caste_SC==1)),
            fixed_effects = ~tehsilid, clusters = villageid, se_type = "stata") 

controls = "top20land + birthyr + vwreserv_gp12 + vscstreserv_gp12"
reds.lm7 = paste(base, "+", controls, "+", "prop_SCnow + prop_STnow") %>% as.formula() %>%
  lm_robust(data = reds.scst.matched %>% filter(period_new==3 & nr_short!=1 & fem==1 & (Caste_ST==1|Caste_SC==1)),
            fixed_effects = ~districtid, clusters = villageid, se_type = "stata")

means = c(
  mean(reds_p3$caste_interact_bi_comp[reds_p3$fem==1 & reds_p3$Caste_ST==1], na.rm = T),
  mean(reds_p3_rand$caste_interact_bi_comp[reds_p3_rand$fem==1 & reds_p3_rand$Caste_ST==1], na.rm = T),
  mean(reds_p3_rand$caste_interact_bi_comp[reds_p3_rand$fem==1 & reds_p3_rand$Caste_ST==1], na.rm = T),
  mean(reds_p3_rand$caste_interact_bi_comp[reds_p3_rand$fem==1 & reds_p3_rand$Caste_ST==1], na.rm = T),
  mean(reds.st.matched$caste_interact_bi_comp[reds.st.matched$period_new==3 & reds.st.matched$nr_short!=1 & reds.st.matched$fem==1 & reds.st.matched$Caste_ST==1], na.rm = T),
  mean(reds_p3_rand$caste_interact_bi_comp[reds_p3_rand$fem==1 & reds_p3_rand$Caste_ST==1], na.rm = T),
  mean(reds.scst.matched$caste_interact_bi_comp[reds.scst.matched$period_new==3 & reds.scst.matched$nr_short!=1 & reds.scst.matched$fem==1 & (reds.scst.matched$Caste_ST==1 | reds.scst.matched$Caste_SC==1)], na.rm = T)
)


screenreg(l = list(reds.lm1, reds.lm2, reds.lm3, 
                   reds.lm4, reds.lm5, reds.lm6, 
                   reds.lm7), include.ci = F, include.rmse = F,
          include.nobs = T,include.adjrs = FALSE, include.rsquared = F,
          digits = 3, include.nclust = F,
          omit.coef = c("Intercept|top20land|birthyr|fem|Caste_ST|Caste_SC|vwreserv_gp12|vstreserv_gp12|prop_STnow|vscstreserv_gp12|prop_SCnow"),
          custom.coef.names = c("Women's Quota", "Caste Quota (ST)", "W X C (ST) Quota", 
                                "Caste Quota (ST/SC)", "W X Caste (ST/SC) Quota"),
          stars = c(0.01, 0.05, 0.1),
          custom.model.names = c("ST Women", "ST Women", "ST Women", "ST Women", "ST Women", "ST Women", "ST Women"),
          custom.gof.rows = list(
            "DV Means" = means,
            "Caste Quota Type" = c("ST", "ST", "ST", "ST", "ST", "ST/SC", "ST/SC"),
            "District FE" = c("No", "Yes", "Yes", "No", "Yes", "No", "Yes"),
            "Subdistrict FE" = c("No", "No", "No", "Yes", "No", "Yes", "No"),
            "Matched Data" = c("No", "No", "No", "No", "Yes", "No", "Yes"),
            "Controls" = c("No", "No", "Yes", "Yes", "Yes", "Yes", "Yes")), 
          file = paste0(tab.out, "tableE9.tex"))


# LONGTERM EFFECTS 
####### TABLE G.16 #############
base = "caste_interact_bi_comp ~ vwreserv_gp12 + vstreserv_gp12 + vwrXvst_gp12"
controls = "top20land + birthyr + Caste_ST + Caste_SC + fem"

reds.lm1 = paste(base) %>% as.formula() %>%
  lm_robust(data = reds_p3, clusters = villageid, se_type = "stata") 

reds.lm2 = paste(base) %>% as.formula() %>%
  lm_robust(data = reds_p3_rand,
            fixed_effects = ~districtid, clusters = villageid, se_type = "stata") 

reds.lm3 = paste(base, "+", controls, "+", "prop_STnow") %>% as.formula() %>%
  lm_robust(data = reds_p3_rand,
            fixed_effects = ~districtid, clusters = villageid, se_type = "stata") 

reds.lm4 = paste(base, "+", controls, "+", "prop_STnow") %>% as.formula() %>%
  lm_robust(data = reds_p3_rand,
            fixed_effects = ~tehsilid, clusters = villageid, se_type = "stata") 

reds.lm5 = paste(base, "+", controls, "+ prop_STnow") %>% as.formula() %>%
  lm_robust(data = reds.st.matched %>% filter(period_new==3 & nr_short!=1),
            fixed_effects = ~districtid,
            clusters = villageid, se_type = "stata") 

base = "caste_interact_bi_comp ~ vwreserv_gp12 + vscstreserv_gp12 + vwrXvscst_gp12"
reds.lm6 = paste(base, "+", controls, "+", "prop_SCnow + prop_STnow") %>% as.formula() %>%
  lm_robust(data = reds_p3_rand,
            fixed_effects = ~tehsilid, clusters = villageid, se_type = "stata") 

controls = "top20land + birthyr + Caste_ST + Caste_SC + fem"
reds.lm7 = paste(base, "+", controls, "+", "prop_STnow + prop_SCnow") %>% as.formula() %>%
  lm_robust(data = reds.scst.matched %>% filter(period_new==3 & nr_short!=1),
            fixed_effects = ~districtid, clusters = villageid, se_type = "stata")



means = c(mean(reds_p3$caste_interact_bi_comp, na.rm = T), 
          mean(reds_p3_rand$caste_interact_bi_comp, na.rm = T), 
          mean(reds_p3_rand$caste_interact_bi_comp, na.rm = T),
          mean(reds_p3_rand$caste_interact_bi_comp, na.rm = T),
          mean(reds.st.matched$caste_interact_bi_comp, na.rm = T),
          mean(reds_p3_rand$caste_interact_bi_comp, na.rm = T),
          mean(reds.scst.matched$caste_interact_bi_comp, na.rm = T))
screenreg(l = list(reds.lm1, reds.lm2, reds.lm3, 
                   reds.lm4, reds.lm5, reds.lm6, 
                   reds.lm7), include.ci = F, include.rmse = F,
          include.nobs = T,include.adjrs = FALSE, include.rsquared = F,
          digits = 3, include.nclust = F,
          omit.coef = c("Intercept|top20land|birthyr|fem|Caste_ST|Caste_SC|prop_STnow|prop_SCnow"),
          custom.coef.names = c("Prior Women's Quota", "Prior Caste Quota (ST)", "Prior W X C (ST) Quota", 
                                "Prior Caste Quota (ST/SC)", "Prior W X Caste (ST/SC) Quota"),
          stars = c(0.01, 0.05, 0.1),
          custom.model.names = c("All", "All", "All", "All", "All", "All", "All"),
          custom.gof.rows = list(
            "Caste Quota Type" = c("ST", "ST", "ST", "ST", "ST", "ST/SC", "ST/SC"),
            "DV Mean" = means,
            "District FE" = c("No", "Yes", "Yes", "No", "Yes", "No", "Yes"),
            "Subdistrict FE" = c("No", "No", "No", "Yes", "No", "Yes", "No"),
            "Matched Data" = c("No", "No", "No", "No", "Yes", "No", "Yes"),
            "Controls" = c("No", "No", "Yes", "Yes", "Yes", "Yes", "Yes")), 
          file = paste0(tab.out, "tableG16.tex"))




### TABLE G.17 #####
base = "caste_interact_bi_comp ~ vwreserv_gp12 + vstreserv_gp12 + vwrXvst_gp12"
controls = "top20land + birthyr + Caste_ST + Caste_SC"

reds.lm1 = paste(base) %>% as.formula() %>%
  lm_robust(data = reds_p3 %>% filter(fem==1), clusters = villageid, se_type = "stata") 

reds.lm2 = paste(base) %>% as.formula() %>%
  lm_robust(data = reds_p3_rand%>% filter(fem==1),
            fixed_effects = ~districtid, clusters = villageid, se_type = "stata") 

reds.lm3 = paste(base, "+", controls, "+", "prop_STnow") %>% as.formula() %>%
  lm_robust(data = reds_p3_rand%>% filter(fem==1),
            fixed_effects = ~districtid, clusters = villageid, se_type = "stata") 

reds.lm4 = paste(base, "+", controls, "+", "prop_STnow") %>% as.formula() %>%
  lm_robust(data = reds_p3_rand%>% filter(fem==1),
            fixed_effects = ~tehsilid, clusters = villageid, se_type = "stata") 

reds.lm5 = paste(base, "+", controls, "+ prop_STnow") %>% as.formula() %>%
  lm_robust(data = reds.st.matched %>% filter(period_new==3 & nr_short!=1 & fem==1),
            fixed_effects = ~districtid,
            clusters = villageid, se_type = "stata") 

base = "caste_interact_bi_comp ~ vwreserv_gp12 + vscstreserv_gp12 + vwrXvscst_gp12"
reds.lm6 = paste(base, "+", controls, "+", "prop_SCnow + prop_STnow") %>% as.formula() %>%
  lm_robust(data = reds_p3_rand%>% filter(fem==1),
            fixed_effects = ~tehsilid, clusters = villageid, se_type = "stata") 

controls = "top20land + birthyr + Caste_ST + Caste_SC"
# reds.scst.matched = reds.scst.matched %>%
#   mutate(caste_interact_bi_comp = ifelse(is.na(caste_interact_bi), 0, caste_interact_bi))
reds.lm7 = paste(base, "+", controls, "+ prop_SCnow + prop_STnow") %>% as.formula() %>%
  lm_robust(data = reds.scst.matched %>% filter(period_new==3 & nr_short!=1 & fem==1),
            fixed_effects = ~districtid, clusters = villageid, se_type = "stata")


means = c(mean(reds_p3$caste_interact_bi_comp[reds_p3$fem==1], na.rm = T), 
          mean(reds_p3_rand$caste_interact_bi_comp[reds_p3_rand$fem==1], na.rm = T), 
          mean(reds_p3_rand$caste_interact_bi_comp[reds_p3_rand$fem==1], na.rm = T),
          mean(reds_p3_rand$caste_interact_bi_comp[reds_p3_rand$fem==1], na.rm = T),
          mean(reds.st.matched$caste_interact_bi_comp[reds.st.matched$fem==1], na.rm = T),
          mean(reds_p3_rand$caste_interact_bi_comp[reds_p3_rand$fem==1], na.rm = T),
          mean(reds.scst.matched$caste_interact_bi_comp[reds.scst.matched$fem==1], na.rm = T))

screenreg(l = list(reds.lm1, reds.lm2, reds.lm3, 
                   reds.lm4, reds.lm5, reds.lm6, 
                   reds.lm7), include.ci = F, include.rmse = F,
          include.nobs = T,include.adjrs = FALSE, include.rsquared = F,
          digits = 3, include.nclust = F,
          omit.coef = c("Intercept|top20land|birthyr|fem|Caste_ST|Caste_SC|prop_STnow|prop_SCnow"),
          custom.coef.names = c("Prior Women's Quota", "Prior Caste Quota (ST)", "Prior W X C (ST) Quota", 
                                "Prior Caste Quota (ST/SC)", "Prior W X Caste (ST/SC) Quota"),
          stars = c(0.01, 0.05, 0.1),
          custom.model.names = c("Women", "Women", "Women", "Women", "Women", "Women", "Women"),
          custom.gof.rows = list(
            "Caste Quota Type" = c("ST", "ST", "ST", "ST", "ST", "ST/SC", "ST/SC"),
            "DV Mean" = means,
            "District FE" = c("No", "Yes", "Yes", "No", "Yes", "No", "Yes"),
            "Subdistrict FE" = c("No", "No", "No", "Yes", "No", "Yes", "No"),
            "Matched Data" = c("No", "No", "No", "No", "Yes", "No", "Yes"),
            "Controls" = c("No", "No", "Yes", "Yes", "Yes", "Yes", "Yes")), 
          file = paste0(tab.out, "tableG17.tex"))




### TABLE G.18 ######
base = "caste_interact_bi_comp ~ vwreserv_gp12 + vstreserv_gp12 + vwrXvst_gp12"
controls = "top20land + birthyr"

reds.lm1 = paste(base) %>% as.formula() %>%
  lm_robust(data = reds_p3 %>% filter(fem==1 & Caste_ST==1), clusters = villageid, se_type = "stata") 

reds.lm2 = paste(base) %>% as.formula() %>%
  lm_robust(data = reds_p3_rand %>% filter(fem==1 & Caste_ST==1),
            fixed_effects = ~districtid, clusters = villageid, se_type = "stata") 

reds.lm3 = paste(base, "+", controls, "+", "prop_STnow") %>% as.formula() %>%
  lm_robust(data = reds_p3_rand %>% filter(fem==1 & Caste_ST==1),
            fixed_effects = ~districtid, clusters = villageid, se_type = "stata") 

reds.lm4 = paste(base, "+", controls, "+", "prop_STnow") %>% as.formula() %>%
  lm_robust(data = reds_p3_rand %>% filter(fem==1 & Caste_ST==1),
            fixed_effects = ~tehsilid, clusters = villageid, se_type = "stata") 

# reds.st.matched = reds.st.matched %>%
#   mutate(caste_interact_bi_comp = ifelse(is.na(caste_interact_bi), 0, caste_interact_bi))
reds.lm5 = paste(base, "+", controls, "+ prop_STnow") %>% as.formula() %>%
  lm_robust(data = reds.st.matched %>% filter(period_new==3 & nr_short!=1 & fem==1 & Caste_ST==1),
            fixed_effects = ~districtid,
            clusters = villageid, se_type = "stata") 

controls = "top20land + birthyr  "
base = "caste_interact_bi_comp ~ vwreserv_gp12 + vscstreserv_gp12 + vwrXvscst_gp12"
reds.lm6 = paste(base, "+", controls, "+", "prop_SCnow + prop_STnow") %>% as.formula() %>%
  lm_robust(data = reds_p3_rand %>% filter(period_new==3 & nr_short!=1 & fem==1 & (Caste_ST==1|Caste_SC==1)),
            fixed_effects = ~tehsilid, clusters = villageid, se_type = "stata") 

controls = "top20land + birthyr"
# reds.scst.matched = reds.scst.matched %>%
#   mutate(caste_interact_bi_comp = ifelse(is.na(caste_interact_bi), 0, caste_interact_bi))
reds.lm7 = paste(base, "+", controls, "+ prop_SCnow + prop_STnow") %>% as.formula() %>%
  lm_robust(data = reds.scst.matched %>% filter(period_new==3 & nr_short!=1 & fem==1 & (Caste_ST==1|Caste_SC==1)),
            fixed_effects = ~districtid, clusters = villageid, se_type = "stata")



means = c(mean(reds_p3$caste_interact_bi_comp[reds_p3$fem==1 & reds_p3$Caste_ST==1], na.rm = T), 
          mean(reds_p3_rand$caste_interact_bi_comp[reds_p3_rand$fem==1 & reds_p3_rand$Caste_ST==1], na.rm = T), 
          mean(reds_p3_rand$caste_interact_bi_comp[reds_p3_rand$fem==1 & reds_p3_rand$Caste_ST==1], na.rm = T),
          mean(reds_p3_rand$caste_interact_bi_comp[reds_p3_rand$fem==1 & reds_p3_rand$Caste_ST==1], na.rm = T),
          mean(reds.st.matched$caste_interact_bi_comp[reds.st.matched$fem==1 & reds.st.matched$Caste_ST==1], na.rm = T),
          mean(reds_p3_rand$caste_interact_bi_comp[reds_p3_rand$fem==1 & (reds_p3_rand$Caste_ST==1|reds_p3_rand$Caste_SC==1)], na.rm = T),
          mean(reds.scst.matched$caste_interact_bi_comp[reds.scst.matched$fem==1 & (reds.scst.matched$Caste_ST==1|reds.scst.matched$Caste_SC==1)], na.rm = T))

screenreg(l = list(reds.lm1, reds.lm2, reds.lm3, 
                   reds.lm4, reds.lm5, reds.lm6, 
                   reds.lm7), include.ci = F, include.rmse = F,
          include.nobs = T,include.adjrs = FALSE, include.rsquared = F,
          digits = 3, include.nclust = F,
          omit.coef = c("Intercept|top20land|birthyr|fem|Caste_ST|Caste_SC|prop_STnow|prop_SCnow"),
          custom.coef.names = c("Prior Women's Quota", "Prior Caste Quota (ST)", "Prior W X C (ST) Quota", 
                                "Prior Caste Quota (ST/SC)", "Prior W X Caste (ST/SC) Quota"),
          stars = c(0.01, 0.05, 0.1),
          custom.model.names = c("ST Women", "ST Women", "ST Women", "ST Women", "ST Women", "ST Women", "ST Women"),
          custom.gof.rows = list(
            "Caste Quota Type" = c("ST", "ST", "ST", "ST", "ST", "ST/SC", "ST/SC"),
            "DV Mean" = means,
            "District FE" = c("No", "Yes", "Yes", "No", "Yes", "No", "Yes"),
            "Subdistrict FE" = c("No", "No", "No", "Yes", "No", "Yes", "No"),
            "Matched Data" = c("No", "No", "No", "No", "Yes", "No", "Yes"),
            "Controls" = c("No", "No", "Yes", "Yes", "Yes", "Yes", "Yes")), 
          file = paste0(tab.out, "tableG18.tex"))


########## TABLE F.10  ##########
# REDS
base = "caste_interact_bi_comp ~ vwreserv*vstreserv"
iv_base = "caste_interact_bi_comp ~ woman_pradhan*ST_pradhan"
controls = "top20land + birthyr + Caste_ST + Caste_SC + fem + vwreserv_gp12 + vstreserv_gp12"

reds.lm1 = paste(iv_base,"+", controls, "+", "prop_STnow", "|", "vwreserv*vstreserv", "+", controls, "+", "prop_STnow") %>% as.formula() %>%
  iv_robust(data = reds_p3_rand %>% filter(nr_short!=1),
            fixed_effects = ~tehsilid, clusters = villageid, se_type = "stata") 
controls = "top20land + birthyr + Caste_ST + Caste_SC + vwreserv_gp12 + vstreserv_gp12"
reds.lm2 = paste(iv_base,"+", controls, "+", "prop_STnow", "|", "vwreserv*vstreserv", "+", controls, "+", "prop_STnow") %>% as.formula() %>%
  iv_robust(data = reds_p3_rand %>% filter(nr_short!=1 & fem==1),
            fixed_effects = ~tehsilid, clusters = villageid, se_type = "stata") 
controls = "top20land + birthyr + vwreserv_gp12 + vstreserv_gp12"

reds.lm3 = paste(iv_base,"+", controls, "+", "prop_STnow", "|", "vwreserv*vstreserv", "+", controls, "+", "prop_STnow") %>% as.formula() %>%
  iv_robust(data = reds_p3_rand%>% filter(nr_short!=1 & fem==1 & Caste_ST==1),
            fixed_effects = ~tehsilid, clusters = villageid, se_type = "stata") 

means = c(mean(reds_p3_rand$caste_interact_bi_comp, na.rm = T), 
          mean(reds_p3_rand$caste_interact_bi_comp[reds_p3_rand$fem==1], na.rm = T), 
          mean(reds_p3_rand$caste_interact_bi_comp[reds_p3_rand$fem==1 & reds_p3_rand$Caste_ST==1], na.rm = T))


screenreg(l = list(reds.lm1, reds.lm2, reds.lm3), include.ci = F, include.rmse = F,
          include.nobs = T,include.adjrs = FALSE, include.rsquared = F,
          digits = 3, include.nclust = F,
          omit.coef = c("Intercept|top20land|birthyr|fem|Caste_ST|Caste_SC|vwreserv_gp12|prop_STnow|vstreserv_gp12|prop_SCnow"),
          custom.coef.names = c("Woman pradhan", "ST pradhan", "W X ST pradhan"),
          stars = c(0.01, 0.05, 0.1),
          custom.model.names = c("All", "Women", "ST Women"),
          custom.gof.rows = list(
            "Caste Quota Type" = c("ST", "ST", "ST"),
            "DV Mean" = means,
            "Subdistrict FE" = c("Yes", "Yes", "Yes"),
            "Matched Data" = c("No", "No", "No"),
            "Controls" = c("Yes", "Yes", "Yes")), 
          file = paste0(tab.out, "tableF10.tex"))





######## TABLE F.11 ######
base = "caste_interact_bi_comp ~ vwreserv*vscstreserv"
iv_base = "caste_interact_bi_comp ~ woman_pradhan*ST_SC_pradhan"
controls = "top20land + birthyr + Caste_ST + Caste_SC + fem+ vwreserv_gp12 + vscstreserv_gp12"

reds.lm4 = paste(iv_base,"+", controls, "+", "prop_STnow", "|", "vwreserv*vscstreserv", "+", controls, "+", "prop_STnow") %>% as.formula() %>%
  iv_robust(data = reds_p3_rand %>% filter(nr_short!=1),
            fixed_effects = ~tehsilid, clusters = villageid, se_type = "stata") 
controls = "top20land + birthyr + Caste_ST + Caste_SC + vwreserv_gp12 + vscstreserv_gp12"

reds.lm5 = paste(iv_base,"+", controls, "+", "prop_STnow", "|", "vwreserv*vscstreserv", "+", controls, "+", "prop_STnow") %>% as.formula() %>%
  iv_robust(data = reds_p3_rand %>% filter(nr_short!=1 & fem==1),
            fixed_effects = ~tehsilid, clusters = villageid, se_type = "stata") 
controls = "top20land + birthyr + vwreserv_gp12 + vscstreserv_gp12"

reds.lm6 = paste(iv_base,"+", controls, "+", "prop_STnow", "|", "vwreserv*vscstreserv", "+", controls, "+", "prop_STnow") %>% as.formula() %>%
  iv_robust(data = reds_p3_rand%>% filter(nr_short!=1 & fem==1 & (Caste_ST==1|Caste_SC==1)),
            fixed_effects = ~tehsilid, clusters = villageid, se_type = "stata") 

means = c(mean(reds_p3_rand$caste_interact_bi_comp, na.rm = T), 
          mean(reds_p3_rand$caste_interact_bi_comp[reds_p3_rand$fem==1], na.rm = T), 
          mean(reds_p3_rand$caste_interact_bi_comp[reds_p3_rand$fem==1 & (reds_p3_rand$Caste_ST==1|reds_p3_rand$Caste_SC==1)], na.rm = T))


screenreg(l = list(reds.lm4, reds.lm5, reds.lm6),
          include.ci = F, include.rmse = F,
          include.nobs = T,include.adjrs = FALSE, include.rsquared = F,
          digits = 3, include.nclust = F,
          omit.coef = c("Intercept|top20land|birthyr|fem|Caste_ST|Caste_SC|vscstreserv_gp12|vwreserv_gp12|prop_STnow|vstreserv_gp12|prop_SCnow"),
          custom.coef.names = c("Woman pradhan", "ST/SC pradhan", "W X ST/SC pradhan"),
          stars = c(0.01, 0.05, 0.1),
          custom.model.names = c("All", "Women", "ST/SC Women"),
          custom.gof.rows = list(
            "Caste Quota Type" = c("ST/SC", "ST/SC", "ST/SC"),
            "DV Mean" = means,
            "Subdistrict FE" = c("Yes", "Yes", "Yes"),
            "Matched Data" = c("No", "No", "No"),
            "Controls" = c("Yes", "Yes", "Yes")), 
          file = paste0(tab.out, "tableF11.tex"))

########## TABLE F.12 ##############

base = "caste_interact_bi ~ vwreserv*vstreserv"
controls = "top20land + birthyr + Caste_SC + Caste_ST + fem + vwreserv_gp12 + vstreserv_gp12"

reds.lm1 = paste(base) %>% as.formula() %>%
  lm_robust(data = reds_p3, clusters = villageid, se_type = "stata") 

reds.lm2 = paste(base) %>% as.formula() %>%
  lm_robust(data = reds_p3_rand,
            fixed_effects = ~districtid, clusters = villageid, se_type = "stata") 

reds.lm3 = paste(base, "+", controls, "+", "prop_STnow") %>% as.formula() %>%
  lm_robust(data = reds_p3_rand,
            fixed_effects = ~districtid, clusters = villageid, se_type = "stata") 

reds.lm4 = paste(base, "+", controls, "+", "prop_STnow") %>% as.formula() %>%
  lm_robust(data = reds_p3_rand,
            fixed_effects = ~tehsilid, clusters = villageid, se_type = "stata") 



reds.lm5 = paste(base, "+", controls) %>% as.formula() %>%
  lm_robust(data = reds.st.matched %>% filter(period_new==3 & nr_short!=1),
            fixed_effects = ~districtid,
            clusters = villageid, se_type = "stata") 

base = "caste_interact_bi ~ vwreserv*vscstreserv"
controls = "top20land + birthyr + Caste_SC + Caste_ST + fem + vwreserv_gp12 + vscstreserv_gp12"

reds.lm6 = paste(base, "+", controls, "+", "prop_SCnow + prop_STnow") %>% as.formula() %>%
  lm_robust(data = reds_p3_rand,
            fixed_effects = ~tehsilid, clusters = villageid, se_type = "stata") 

controls = "top20land + birthyr + Caste_SC + Caste_ST + fem + vwreserv_gp12 + vscstreserv_gp12"

reds.lm7 = paste(base, "+", controls, "+", "prop_SCnow + prop_STnow") %>% as.formula() %>%
  lm_robust(data = reds.scst.matched %>% filter(period_new==3 & nr_short!=1),
            fixed_effects = ~districtid, clusters = villageid, se_type = "stata")


means = c(mean(reds_p3$caste_interact_bi, na.rm = T), 
          mean(reds_p3_rand$caste_interact_bi, na.rm = T), 
          mean(reds_p3_rand$caste_interact_bi, na.rm = T),
          mean(reds_p3_rand$caste_interact_bi, na.rm = T),
          mean(reds.st.matched$caste_interact_bi, na.rm = T),
          mean(reds_p3_rand$caste_interact_bi, na.rm = T),
          mean(reds.scst.matched$caste_interact_bi, na.rm = T))


screenreg(l = list(reds.lm1, reds.lm2, reds.lm3, reds.lm4, reds.lm5, 
                   reds.lm6, reds.lm7),
          include.ci = F, include.rmse = F,
          include.nobs = T,include.adjrs = FALSE, include.rsquared = F,
          digits = 3, include.nclust = F,
          omit.coef = c("Intercept|top20land|birthyr|fem|Caste_ST|Caste_SC|vscstreserv_gp12|vwreserv_gp12|prop_STnow|vstreserv_gp12|prop_SCnow"),
          custom.coef.names = c("Women's Quota", "Caste (ST) Quota", "W X C (ST) Quota", 
                                "Caste (ST/SC) Quota", "W X C (ST/SC) Quota"),
          stars = c(0.01, 0.05, 0.1),
          custom.model.names = c("All", "All", "All", "All", "All", "All", "All"),
          custom.gof.rows = list(
            "Caste Quota Type" = c("ST", "ST", "ST", "ST", "ST", "ST/SC", "ST/SC"),
            "DV Mean" = means,
            "District FE" = c("No", "Yes", "Yes", "No", "Yes", "No", "Yes"),
            "Sub-district FE" = c("No", "No", "No", "Yes", "No", "Yes", "No"),
            "Matched Data" = c("No", "No", "No", "No", "Yes", "No", "Yes"),
            "Controls" = c("No", "No", "Yes", "Yes", "Yes", "Yes", "Yes")), 
          file = paste0(tab.out, "tableF12.tex"))


############# TABLE F.13 ########################
# make sure outcome is coded as a factor
reds_p3_rand$caste_interact_sc_comp_fact = factor(reds_p3_rand$caste_interact_sc_comp, ordered = T, levels = 1:5)

controls = "top20land + birthyr + fem + Caste_ST + Caste_SC + vwreserv_gp12 + vstreserv_gp12"
base = paste0("caste_interact_sc_comp_fact ~ vwreserv*vstreserv + prop_STnow")

m6 = polr(base %>% as.formula, data = reds_p3_rand, Hess = T)
m6_coef = exp(coefficients(m6))
m6_t = coeftest(m6)[,3]


controls = "top20land + birthyr + Caste_ST + Caste_SC + vwreserv_gp12 + vstreserv_gp12"
m8 = polr(base %>% as.formula, data = reds_p3_rand %>%
            dplyr::filter(fem==1), Hess = T)
m8_coef = exp(coefficients(m8))
m8_t = coeftest(m8)[,3]


controls = "top20land + birthyr + vwreserv_gp12 + vstreserv_gp12"

m10 = polr(base %>% as.formula, data = reds_p3_rand %>%
             dplyr::filter(Caste_ST==1 & fem==1), Hess = T)
m10_coef = exp(coefficients(m10))
m10_t = coeftest(m10)[,3]


stargazer(m6, m8, m10,  report = c("vc*t"), star.cutoffs = NA,
          coef = list( m6_coef, m8_coef, m10_coef),
          booktabs = T, t = list(m6_t, m8_t, m10_t),
          covariate.labels = c("Women's Quota", "Caste Quota", "W X C Quota"),
          title = "Impact of Reservations on Inter-Group Relations, Multinomial Logit Model",
          omit = c("top20land|educghi|educmghi|cohort_young|cohort_mid|cohort_old|prop_STnow|Intercept|vwreserv_gp12|vstreserv_gp12"),
          add.lines = list(
            c("Caste Quota Type", "ST", "ST", "ST"),
            c("Model", "Ordinal logit", "Ordinal logit", "Ordinal logit"),
            c("Sample", "Entire", "All Women", "ST Women"),
            c("District FE", "No", "No", "No"), 
            c("Controls", "Yes", "Yes", "Yes")
          ),
          out = paste0(tab.out, "tableF13.tex"))


########### TABLE F.14 ##########
reds.lmer1 = lmer(caste_interact_bi_comp ~ vwreserv*vstreserv + top20land + birthyr + fem + Caste_ST + Caste_SC +
                    vwreserv_gp12 + vwreserv_gp12 + 
                    prop_STnow + (1 | tehsilid) + (1 | districtid), data = reds_p3)

reds.lmer2 = lmer(caste_interact_bi_comp ~ vwreserv*vstreserv + top20land + birthyr + fem + Caste_ST + Caste_SC +
                    vwreserv_gp12 + vwreserv_gp12 + 
                    prop_STnow + (1 | tehsilid) + (1 | districtid), data = reds_p3_rand)

reds.lmer3 = lmer(caste_interact_bi_comp ~ vwreserv*vstreserv + top20land + birthyr + fem + Caste_ST + Caste_SC +
                    vwreserv_gp12 + vwreserv_gp12 + 
                    prop_STnow + (1 | villageid) + (1 | tehsilid) + (1 | districtid), data = reds_p3_rand)

reds.lmer4 = lmer(caste_interact_bi_comp ~ vwreserv*vstreserv + top20land + birthyr + Caste_ST + Caste_SC +
                    vwreserv_gp12 + vwreserv_gp12 + 
                    prop_STnow + (1 | villageid) + (1 | tehsilid) + (1 | districtid), data = reds_p3_rand %>% filter(fem==1))

reds.lmer5 = lmer(caste_interact_bi_comp ~ vwreserv*vstreserv + top20land + birthyr +
                    vwreserv_gp12 + vwreserv_gp12 + 
                    prop_STnow + (1 | villageid) + (1 | tehsilid) + (1 | districtid), data = reds_p3_rand %>% filter(fem==1 & Caste_ST==1))

reds.lmer6 = lmer(caste_interact_bi_comp ~ vwreserv*vscstreserv + top20land + birthyr +
                    vwreserv_gp12 + vwreserv_gp12 + 
                    prop_STnow + (1 | villageid) + (1 | tehsilid) + (1 | districtid), data = reds_p3_rand %>% filter(fem==1 & (Caste_ST==1|Caste_SC==1)))


screenreg(l = list(reds.lmer1, reds.lmer2, reds.lmer3, reds.lmer4, reds.lmer5, reds.lmer6),
               caption = "\\textbf{Effect of Reservations on Attitude Towards Other Castes, All-India, ST Women, Multilevel Models}",
               caption.above = T, digits = 3, booktabs = T,
               naive = T, include.ci = F,
               omit.coef = c("age|vote_panch|fem|birthyr|Caste_ST|Caste_SC|tehsilid|districtid|period_new|top20land|educghi|educmghi|cohort_young|cohort_mid|cohort_old|Intercept|prop_SCnow|prop_STnow|prop_st|vwreserv_gp12|vstreserv_gp12"),
               custom.header = list("Interacting With Other Castes" = 1:6),
               custom.coef.names = c("Women's Quotas", "ST Quotas", "W. $\\times$ ST Q.", "SC/ST Quotas", "W $\\times$ SC/ST Quotas"),
               custom.gof.rows = list("Sample" = c("All", "All", "All", "All Women", "ST Women", "ST/SC Women"),
                                      "States" = c("All", "Randomly.Assigned", "RA", "RA", "RA","RA"),
                                      "Controls" = c("Yes","Yes","Yes", "Yes", "Yes", "Yes")),
               stars = c(0.01, 0.05, 0.1), 
               custom.note = paste0("%stars"), file = paste0(tab.out, "tableF14.tex"))




######### TABLE F.15  ###########
base = "caste_interact_bi_comp ~ vwreserv*vstreserv"
controls = "top20land + birthyr + fem + Caste_ST + Caste_SC + vwreserv_gp12 + vstreserv_gp12"

reds.lm4 = paste(base, "+", controls, "+", "prop_STnow ", "+ as.factor(tehsilid)") %>%
  as.formula() %>%
  lm(data = reds_p3_rand) 
# ST proportion - W X ST Quota
quota.sensitivity <- sensemakr(model = reds.lm4, 
                                treatment = "vwreserv:vstreserv",
                                benchmark_covariates = "prop_STnow",
                                kd = 1:3,
                                ky = 1:3, 
                                q = 1,
                                alpha = 0.05, 
                                reduce = TRUE)
tab <- summary(quota.sensitivity) %>%
  data.frame() %>%
  dplyr::select(-R2dz.x, -R2yz.dx) %>%
  xtable(digits = 4)
print(tab,file = paste0(tab.out, "tableF15_panel1.tex"))
# Female - W X ST Quota
quota.sensitivity.fem <- sensemakr(model = reds.lm4, 
                               treatment = "vwreserv:vstreserv",
                               benchmark_covariates = "fem",
                               kd = 1:3,
                               ky = 1:3, 
                               q = 1,
                               alpha = 0.05, 
                               reduce = TRUE)


tab <- summary(quota.sensitivity.fem)  %>%
  data.frame() %>%
  dplyr::select(-R2dz.x, -R2yz.dx) %>%
  xtable(digits = 4)
print(tab, file = paste0(tab.out, "tableF15_panel2.tex"))

